function WS = willmott_score(modeldata, obsdata)
%
% this computes the willmott score based on two time series
%
% WS = 1 - MSE / <(|m - <o>| + |o - <o>|)^2> 
%  
% where MSE is mean squared error = <(m - o)>^2
% and m = modeled, o = observed,
% and || is abs. value, and <> is mean
%************************************************

m = modeldata;
o = obsdata;

if(length(m)~=length(o))
    error('Vectors must be the same length!');
end

MSE = nanmean((m - o).^2);
denom1 = abs(m - nanmean(o));
denom2 = abs(o - nanmean(o));
denom = nanmean((denom1 + denom2).^2);

WS = 1 - MSE./denom;




